Generalized linear mixed effects models were appropriate for this
analysis because the data contain sources of non-independence.
Individual plants were located in the same site as other individuals and
were likely genetically similar. Therefore, we considered site a random
effect.
The response variable is total seeds, which is discrete and non-zero.
The interaction between deviation and the early/late factor was a fixed
effect, as well as the interaction between plot treatment and
conspecific density. Flower number is included (but not assessed) in the
interaction terms to account for the effect of the number of flowers on
number of seeds. We would have included soil moisture as a covariate and
fixed effect in our model, but preliminary analyses involving model
selection did not identify soil moisture as a significant driver of seed
production. Each species is analyzed separately.
Mertensia
We started by limiting the individuals to the open treatment
(excluding the hand-pollinated treatment) and by checking the
distribution of the total seed counts.
mert.open<-mert[(mert$treat=="open"),] #removed hand pollination treatment
plot(mert.nbinom) #plotting to see the fit
The Mertensia seed data fits a negative binomial distribution. The
glmmTMB package is a good package for modeling data with a negative
binomial. The glmmTMB package is also suitable for data that could be
overdispersed or zero-inflated.
Main Model
#mert.glmmtmb<-glmmTMB(dev.seed ~ flowers:(deviation*early.late) + flowers:number.conspecifics/plot.treat + (1|site), family = "nbinom2", data = mert.open) #putting count data into glmmTMB for DHARMa package
#summary(mert.glmmtmb) #summary of glmmTMB model
#mert.counts.simulation<-simulateResiduals(fittedModel = mert.glmmtmb) #creating simulation for Mertensia mixed effects model with DHARMa package
#plot(mert.counts.simulation) #plotting simulation
#testDispersion(mert.counts.simulation) #checking for overdispersion
#testZeroInflation(mert.counts.simulation) #checking for zero inflation
We detected a significant effect of early vs. late blooming on the
total number of seeds in Mertensia individuals, but when we checked for
overdispersion and zero-inflation in the data using the DHARMa package,
the data were slightly zero-inflated. To address this, we adjusted our
glmmtmb model to include zero-inflation.
mert.glmmtmb<-glmmTMB(dev.seed ~ flowers:(deviation*early.late) + flowers:number.conspecifics/plot.treat + (1|site), ziformula=~1, family = "nbinom2", data = mert.open) #putting count data into glmmTMB for DHARMa package
summary(mert.glmmtmb) #summary of glmmTMB model
## Family: nbinom2 ( log )
## Formula:
## dev.seed ~ flowers:(deviation * early.late) + flowers:number.conspecifics/plot.treat +
## (1 | site)
## Zero inflation: ~1
## Data: mert.open
##
## AIC BIC logLik deviance df.resid
## NA NA NA NA 38
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.06364 0.2523
## Number of obs: 47, groups: site, 4
##
## Overdispersion parameter for nbinom2 family (): 2.91
##
## Conditional model:
## Estimate Std. Error z value
## (Intercept) 2.400e+00 2.417e-01 9.929
## flowers:deviation 2.565e-03 1.523e-03 1.684
## flowers:early.latelate 4.751e-02 1.922e-02 2.471
## flowers:number.conspecifics -7.164e-05 8.747e-06 -8.191
## flowers:deviation:early.latelate -1.084e-02 7.623e-03 -1.423
## flowers:number.conspecifics:plot.treatmanip 2.311e-04 3.138e-04 0.736
## Pr(>|z|)
## (Intercept) < 2e-16 ***
## flowers:deviation 0.0921 .
## flowers:early.latelate 0.0135 *
## flowers:number.conspecifics 2.59e-16 ***
## flowers:deviation:early.latelate 0.1548
## flowers:number.conspecifics:plot.treatmanip 0.4615
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.4634 0.5713 -4.312 1.62e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mert.counts.simulation<-simulateResiduals(fittedModel = mert.glmmtmb) #creating simulation for Mertensia mixed effects model with DHARMa package
plot(mert.counts.simulation) #plotting simulation
testDispersion(mert.counts.simulation) #checking for overdispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## ratioObsSim = 0.82711, p-value = 0.44
## alternative hypothesis: two.sided
testZeroInflation(mert.counts.simulation) #checking for zero inflation
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.0593, p-value = 1
## alternative hypothesis: two.sided
When number of flowers is taken into account, seed set for Mertensia
individuals may be affected by blooming before or after the population
peak. Zero-inflation is accounted for in the new model, and the model
for Mertensia count data is no longer zero-inflated.
To figure out the direction of the effect, we looked at the plots
(see “Figures”). Blooming after the population peak may increase
seed set for Mertensia individuals.
Pollen Limitation
Model
To test the pollen limitation individuals, we created a separate
fixed effects model. This model has site and plot treatment as random
effects and the interaction between pollen limitation treatment and
number of flowers as a fixed effect. We included this interaction
because the difference in seed set between the open and hand treatments
is heavily dependent on the number of flowers. The interaction accounts
for the wide variation of number of flowers and ultimately developed
seeds. The response variable is still total number of developed seeds in
Mertensia individuals.
#mert.treat.glmmtmb<-glmmTMB(dev.seed ~ treat/flowers + (1|site/plot.treat), family = "nbinom2", data = mert) #Mertensia model for negative binomial data with glmmTMB package
#summary(mert.treat.glmmtmb) #summary of model output
#mert.counts.treat.simulation<-simulateResiduals(fittedModel = mert.treat.glmmtmb) #creating simulation for Mertensia mixed effects model with DHARMa package
#plot(mert.counts.treat.simulation) #plotting simulation
#testDispersion(mert.counts.treat.simulation) #checking for overdispersion
#testZeroInflation(mert.counts.treat.simulation) #checking for zero inflation
We detected an effect of pollen limitation treatment on Mertensia
total developed seeds, but we had to test the model assumptions of our
pollen limitation model as well. The data in the pollen limitation model
for Mertensia are not overdispersed, but they are zero-inflated. We then
adjusted our model to include zero-inflation and ran the overdispersion
and zero-inflation tests again.
mert.treat.glmmtmb<-glmmTMB(dev.seed ~ treat/flowers + (1|site/plot.treat), ziformula=~1, family = "nbinom2", data = mert) #Mertensia model for negative binomial data with glmmTMB package
summary(mert.treat.glmmtmb) #summary of model output
## Family: nbinom2 ( log )
## Formula: dev.seed ~ treat/flowers + (1 | site/plot.treat)
## Zero inflation: ~1
## Data: mert
##
## AIC BIC logLik deviance df.resid
## 807.4 828.2 -395.7 791.4 92
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## plot.treat:site (Intercept) 0.02086 0.1444
## site (Intercept) 0.03943 0.1986
## Number of obs: 100, groups: plot.treat:site, 8; site, 4
##
## Overdispersion parameter for nbinom2 family (): 4.51
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.255095 0.225157 10.016 < 2e-16 ***
## treatopen-hand -0.068716 0.268730 -0.256 0.798
## treatopen:flowers 0.024104 0.004384 5.498 3.83e-08 ***
## treatopen-hand:flowers 0.027786 0.004423 6.282 3.35e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.2355 0.5407 -5.984 2.18e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mert.counts.treat.simulation<-simulateResiduals(fittedModel = mert.treat.glmmtmb) #creating simulation for Mertensia mixed effects model with DHARMa package
plot(mert.counts.treat.simulation) #plotting simulation
testDispersion(mert.counts.treat.simulation) #checking for overdispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## ratioObsSim = 0.80744, p-value = 0.232
## alternative hypothesis: two.sided
testZeroInflation(mert.counts.treat.simulation) #checking for zero inflation
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 0.97752, p-value = 1
## alternative hypothesis: two.sided
The data are no longer zero-inflated, and we still detected an effect
of pollen limitation on the Mertensia seed counts. Mertensia
individuals are pollen-limited.
Figures
Below is the plot of seed set vs. relative position of blooming in
the population for Mertensia individuals. Plot treatment is indicated by
individual color.
ggplot(mert.open, aes(x=mert.open$relative.position, y=mert.open$dev.seed, group=mert.open$plot.treat)) +
labs(y="Total Developed Seeds",x="Days Relative to Population Peak Bloom", title="Effect of Mertensia Blooming Time on Total Developed Seeds") +
geom_point(aes(color=plot.treat),size=3, alpha = 0.8) +
scale_color_brewer(palette="Dark2") + labs(color="Plot Treatment") +
scale_color_manual(labels= c("control","accelerated"), values = c("control"="darkred","manip"="darkorange")) +
geom_smooth()
ggplot(mert.open, aes(x=mert.open$relative.position, y=mert.open$dev.seed, group=mert.open$early.late)) +
labs(y="Total Developed Seeds",x="Days Relative to Population Peak Bloom", title="Effect of Mertensia Blooming Time on Total Developed Seeds") +
theme_classic() +
geom_point(aes(color=early.late),size=3, alpha = 0.8) +
scale_color_brewer(palette="Dark2") + labs(color="Before or After Peak") +
scale_color_manual(labels= c("before peak","after peak"), values = c("early"="darkred","late"="darkorange")) +
geom_smooth(method = "glm", aes(color=early.late))
We then plotted the interaction between deviation and early/late
instead of relative position in order to fit a linear model.
ggplot(mert.open, aes(x=mert.open$deviation, y=mert.open$dev.seed/mert.open$flowers, group= mert.open$early.late)) +
labs(y="Total Developed Seeds per Flower",x="Days Relative to Population Peak Bloom", title="Effect of Mertensia Deviation from Peak on Total Developed Seeds") +
theme_classic() +
scale_x_continuous(breaks=seq(0,6,1)) +
geom_point(size=3, alpha = 0.8, aes(color=early.late)) +
geom_smooth(method="glm",aes(color=early.late)) +
scale_color_brewer(palette="Dark2") +
labs(color="Before or After Peak") +
scale_color_manual(labels= c("before peak","after peak"), values = c("early"="darkred","late"="darkorange"))
Then we plotted the effect of conspecific density and plot treatment
on total developed seeds.
ggplot(mert.open, aes(x=mert.open$number.conspecifics, y=mert.open$dev.seed/mert.open$flowers, group= mert.open$plot.treat)) +
labs(y="Total Developed Seeds per Flower",x="Conspecific Density", title="Effect of Mertensia Conspecific Density on Total Developed Seeds") +
theme_classic() +
geom_point(size=3, alpha = 0.8, aes(color=plot.treat)) +
geom_smooth(method="glm",aes(color=plot.treat)) +
scale_color_brewer(palette="Dark2") +
labs(color="Plot Treatment") +
scale_color_manual(labels= c("control","accelerated snowmelt"), values = c("control"="darkgreen","manip"="blue"))
Finally, we plotted all of the data with the pollinator active
period. The pollinator active period was the time during which we saw
interactions between pollinators and our species. The active period is
shaded.
mert$midpoint2<-format(as.Date(mert$midpoint), format="%m-%d")
ggplot(mert, aes(x=mert$midpoint2, y=mert$dev.seed/mert$flowers, group= mert$treat)) +
labs(y="Total Developed Seeds per Flower",x="Date", title="Effect of Mertensia Blooming Time and Pollinator Activity on Total Developed Seeds") +
theme_classic() +
geom_point(size=3, alpha = 0.8, aes(color=treat)) +
geom_smooth(aes(color=treat)) +
scale_color_brewer(palette="Dark2") +
labs(color="Open or Hand-pollinated") +
scale_color_manual(labels= c("open","open-hand"), values = c("open"="darkred","open-hand"="darkorange")) +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
annotate("rect",xmin="06-19", xmax="07-03",ymin=0,ymax=Inf,alpha=0.1)
Delphinium
For Delphinium, we again limited the individuals to open treatment
individuals.
delph.open<-delph[(delph$treat=="open"),] #removed hand pollination treatment
plot(delph.nbinom) #plotting distribution fit
The Delphinium data fits a negative binomial distribution. We used
the glmmTMB package again.
Main Model
delph.glmmtmb<-glmmTMB(total.seeds ~ X.flowers:(deviation*early.late) + X.flowers:number.conspecifics/plot.treat + (1|site), family = "nbinom2", data = delph.open) #Delphinium model for negative binomial data with glmmTMB package
#removed treat fixed effect (only open)
summary(delph.glmmtmb) #summary of model output
## Family: nbinom2 ( log )
## Formula:
## total.seeds ~ X.flowers:(deviation * early.late) + X.flowers:number.conspecifics/plot.treat +
## (1 | site)
## Data: delph.open
##
## AIC BIC logLik deviance df.resid
## 403.9 417.6 -193.9 387.9 33
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.02556 0.1599
## Number of obs: 41, groups: site, 4
##
## Overdispersion parameter for nbinom2 family (): 1.37
##
## Conditional model:
## Estimate Std. Error z value
## (Intercept) 2.990040 0.369707 8.088
## X.flowers:deviation 0.026479 0.011771 2.249
## X.flowers:early.latelate 0.137212 0.168217 0.816
## X.flowers:number.conspecifics 0.001297 0.001946 0.666
## X.flowers:deviation:early.latelate -0.010662 0.036032 -0.296
## X.flowers:number.conspecifics:plot.treatmanip -0.002071 0.001801 -1.150
## Pr(>|z|)
## (Intercept) 6.09e-16 ***
## X.flowers:deviation 0.0245 *
## X.flowers:early.latelate 0.4147
## X.flowers:number.conspecifics 0.5051
## X.flowers:deviation:early.latelate 0.7673
## X.flowers:number.conspecifics:plot.treatmanip 0.2501
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
delph.counts.simulation<-simulateResiduals(fittedModel = delph.glmmtmb) #creating simulation for Delphinium mixed effects model with DHARMa package
plot(delph.counts.simulation) #plotting simulation
testDispersion(delph.counts.simulation) #checking for overdispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## ratioObsSim = 0.79078, p-value = 0.448
## alternative hypothesis: two.sided
testZeroInflation(delph.counts.simulation) #checking for zero inflation
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 4.2735, p-value = 0.136
## alternative hypothesis: two.sided
When flower number is taken into account, deviation from the
population peak has a significant effect on Delphinium developed seed
count. We tested the model assumptions as we did with Mertensia seed
count data. The model for the Delphinium count data is not overdispersed
or zero-inflated.
To determine the direction of deviation effect, we looked at the
plot. When flower number is taken into account, seed set of
Delphinium individuals increases with deviation from the population
peak.
Pollen Limitation
Model
Again, we tested the pollen limitation in a separate mixed effects
model.
#delph.treat.glmmtmb<-glmmTMB(total.seeds ~ treat/X.flowers + (1|site/plot.treat), family = "nbinom2", data = delph) #Mertensia model for negative binomial data with glmmTMB package
#summary(delph.treat.glmmtmb) #summary of model output
#delph.counts.treat.simulation<-simulateResiduals(fittedModel = delph.treat.glmmtmb) #creating simulation for Delphinium mixed effects model with DHARMa package
#plot(delph.counts.treat.simulation) #plotting simulation
#testDispersion(delph.counts.treat.simulation) #checking for overdispersion
#testZeroInflation(delph.counts.treat.simulation) #checking for zero inflation
We detected an effect of pollen limitation on the total developed
seeds of Delphinium individuals. We again checked for overdispersion and
zero-inflation in the pollen limitation data. The Delphinum pollen
limitation model is zero-inflated. We adjusted our model to include the
zero-inflation formula.
delph.treat.glmmtmb<-glmmTMB(total.seeds ~ treat/X.flowers + (1|site/plot.treat), ziformula=~1, family = "nbinom2", data = delph) #Mertensia model for negative binomial data with glmmTMB package
summary(delph.treat.glmmtmb) #summary of model output
## Family: nbinom2 ( log )
## Formula: total.seeds ~ treat/X.flowers + (1 | site/plot.treat)
## Zero inflation: ~1
## Data: delph
##
## AIC BIC logLik deviance df.resid
## 811.2 830.6 -397.6 795.2 76
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## plot.treat:site (Intercept) 2.415e-09 4.914e-05
## site (Intercept) 6.478e-02 2.545e-01
## Number of obs: 84, groups: plot.treat:site, 7; site, 4
##
## Overdispersion parameter for nbinom2 family (): 1.97
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.51490 0.40312 6.239 4.42e-10 ***
## treatopen-hand 0.69778 0.50148 1.391 0.164096
## treatopen:X.flowers 0.27607 0.07661 3.603 0.000314 ***
## treatopen-hand:X.flowers 0.16104 0.06725 2.395 0.016636 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.1192 0.5794 -5.383 7.31e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
delph.counts.treat.simulation<-simulateResiduals(fittedModel = delph.treat.glmmtmb) #creating simulation for Delphinium mixed effects model with DHARMa package
plot(delph.counts.treat.simulation) #plotting simulation
testDispersion(delph.counts.treat.simulation) #checking for overdispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## ratioObsSim = 0.79907, p-value = 0.352
## alternative hypothesis: two.sided
testZeroInflation(delph.counts.treat.simulation) #checking for zero inflation
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.065, p-value = 1
## alternative hypothesis: two.sided
The new model is not zero-inflated. We still see an effect of pollen
limitation on the Delphinium count data. Delphinium individuals
are pollen-limited.
Figures
Below are the plots for the Delphinium individuals.
ggplot(delph.open, aes(x=delph.open$relative.position, y=delph.open$total.seeds, group=plot.treat)) +
labs(title="Effect of Delphinium Blooming Time on Total Developed Seeds", y="Total Developed Seeds",x="Days Relative to Population Peak Bloom") +
geom_point(aes(color=plot.treat),size=3, alpha = 0.8) +
scale_color_brewer(palette="Dark2") + labs(color="Plot Treatment") +
scale_color_manual(labels= c("control","accelerated"), values = c("control"="darkred","manip"="darkorange")) +
geom_smooth()
ggplot(delph.open, aes(x=delph.open$relative.position, y=delph.open$total.seeds, group=delph.open$early.late)) +
labs(y="Total Developed Seeds",x="Days Relative to Population Peak Bloom", title="Effect of Delphinium Blooming Time on Total Developed Seeds") +
theme_classic() +
geom_point(aes(color=early.late),size=3, alpha = 0.8) +
scale_color_brewer(palette="Dark2") + labs(color="Before or After Peak") +
scale_color_manual(labels= c("before peak","after peak"), values = c("early"="darkred","late"="darkorange")) +
geom_smooth(method = "glm", aes(color=early.late))
ggplot(delph.open, aes(x=delph.open$deviation, y=(delph.open$total.seeds/delph.open$X.flowers), group= delph.open$early.late)) +
labs(y="Total Developed Seeds per Flower",x="Days Relative to Population Peak Bloom", title="Effect of Delphinium Deviation from Peak on Total Developed Seeds") +
theme_classic() +
scale_x_continuous(breaks=seq(0,10,1)) +
geom_point(size=3, alpha = 0.8, aes(color=early.late)) +
geom_smooth(method="glm",aes(color=early.late)) +
scale_color_brewer(palette="Dark2") +
labs(color="Before or After Peak") +
scale_color_manual(labels= c("before peak","after peak"), values = c("early"="darkred","late"="darkorange"))
ggplot(delph.open, aes(x=delph.open$number.conspecifics, y=delph.open$total.seeds/delph.open$X.flowers, group= delph.open$plot.treat)) +
labs(y="Total Developed Seeds per Flower",x="Conspecific Density", title="Effect of Delphinium Conspecific Density on Total Developed Seeds") +
theme_classic() +
geom_point(size=3, alpha = 0.8, aes(color=plot.treat)) +
geom_smooth(method="glm",aes(color=plot.treat)) +
scale_color_brewer(palette="Dark2") +
labs(color="Plot Treatment") +
scale_color_manual(labels= c("control","accelerated snowmelt"), values = c("control"="darkgreen","manip"="blue"))
Finally, we plotted the pollinator activity period with the total
number of developed seeds and midpoint of individuals.
delph$midpoint2<-format(as.Date(delph$midpoint), format="%m-%d")
ggplot(delph, aes(x=delph$midpoint2, y=delph$total.seeds/delph$X.flowers, group= delph$treat)) +
labs(y="Total Developed Seeds per Flower",x="Date", title="Effect of Delphinium Blooming and Pollinator Activity on Total Developed Seeds") +
theme_classic() +
geom_point(size=3, alpha = 0.8, aes(color=treat)) +
geom_smooth(aes(color=treat)) +
scale_color_brewer(palette="Dark2") +
labs(color="Open or Hand-pollinated") +
scale_color_manual(labels= c("open","open-hand"), values = c("open"="darkred","open-hand"="darkorange")) +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
annotate("rect",xmin="06-27", xmax="07-18",ymin=0,ymax=Inf,alpha=0.1)
Potentilla
We began by limiting the data to the open treatment Potentilla
individuals.
pot.open<-pot[(pot$treat=="open"),] #removed hand pollination treatment
plot(pot.nbinom) #plotting distribution fit
The Potentilla seed data fits a negative binomial distribution. We used
the glmmTMB package again.
Main Model
#pot.glmmtmb<-glmmTMB(total.seeds ~ X.flowers:(deviation * early.late) + X.flowers:number.conspecifics/plot.treat + (1|site), family = "nbinom2", data = pot.open) #Potentilla model for negative binomial data with glmmTMB package
#summary(pot.glmmtmb) #summary of model output
#pot.counts.simulation<-simulateResiduals(fittedModel = pot.glmmtmb) #creating simulation for Potentilla mixed effects model with DHARMa package
#plot(pot.counts.simulation) #plotting simulation
#testDispersion(pot.counts.simulation) #checking for overdispersion
#testZeroInflation(pot.counts.simulation) #checking for zero inflation
The Potentilla data are not overdispersed, but they are
zero-inflated. To fix this, we had to include zero inflation in our
model and again check for zero-inflation.
pot.glmmtmb<-glmmTMB(total.seeds ~ X.flowers:(deviation * early.late) + X.flowers:number.conspecifics/plot.treat + (1|site), ziformula=~1, family = "nbinom2", data = pot.open) #new Potentilla model for negative binomial data with glmmTMB package
summary(pot.glmmtmb) #summary of model output
## Family: nbinom2 ( log )
## Formula:
## total.seeds ~ X.flowers:(deviation * early.late) + X.flowers:number.conspecifics/plot.treat +
## (1 | site)
## Zero inflation: ~1
## Data: pot.open
##
## AIC BIC logLik deviance df.resid
## 280.7 291.4 -131.4 262.7 15
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## site (Intercept) 3.537e-10 1.881e-05
## Number of obs: 24, groups: site, 4
##
## Overdispersion parameter for nbinom2 family (): 2.95
##
## Conditional model:
## Estimate Std. Error z value
## (Intercept) 3.657e+00 4.095e-01 8.930
## X.flowers:deviation 1.220e-02 4.421e-03 2.760
## X.flowers:early.latelate 7.605e-02 2.733e-02 2.783
## X.flowers:number.conspecifics 1.631e-05 3.236e-04 0.050
## X.flowers:deviation:early.latelate -1.423e-02 4.629e-03 -3.074
## X.flowers:number.conspecifics:plot.treatmanip 4.595e-04 3.227e-04 1.424
## Pr(>|z|)
## (Intercept) < 2e-16 ***
## X.flowers:deviation 0.00578 **
## X.flowers:early.latelate 0.00539 **
## X.flowers:number.conspecifics 0.95981
## X.flowers:deviation:early.latelate 0.00211 **
## X.flowers:number.conspecifics:plot.treatmanip 0.15444
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.138 1.024 -3.065 0.00218 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
pot.counts.simulation<-simulateResiduals(fittedModel = pot.glmmtmb) #creating simulation for Potentilla mixed effects model with DHARMa package
plot(pot.counts.simulation) #plotting simulation
testDispersion(pot.counts.simulation) #checking for overdispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## ratioObsSim = 0.8047, p-value = 0.544
## alternative hypothesis: two.sided
testZeroInflation(pot.counts.simulation) #checking for zero inflation
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.0776, p-value = 1
## alternative hypothesis: two.sided
When we incorporated zero-inflation into the model, we found a
significant effect of deviation from the population peak, blooming
before vs after the peak, and the interaction between the two. The
direction of the effect, however, must be determined from the plots (see
“Figures”).
After looking at the plots, we see that when we account for
the number of flowers, blooming after the peak lowers seed set and
deviation has a more negative effect on the late blooming Potentilla
individuals.
Pollen Limitation
Model
We assessed pollen limitation for Potentilla in a mixed effects
model.
#pot.treat.glmmtmb<-glmmTMB(total.seeds ~ treat/X.flowers + (1|site/plot.treat), family = "nbinom2", data = pot) #Potentilla model for negative binomial data with glmmadmb package
#summary(pot.treat.glmmtmb) #summary of model output
#pot.counts.treat.simulation<-simulateResiduals(fittedModel = pot.treat.glmmtmb) #creating simulation for Potentilla mixed effects model with DHARMa package
#plot(pot.counts.treat.simulation) #plotting simulation
#testDispersion(pot.counts.treat.simulation) #checking for overdispersion
#testZeroInflation(pot.counts.treat.simulation) #checking for zero inflation
We detected a significant effect of the pollen limitation treatment
on the total number of developed seeds for Potentilla individuals, but
the Potentilla pollen limitation data were zero-inflated. We again
included the zero-inflation formula into our models.
pot.treat.glmmtmb<-glmmTMB(total.seeds ~ treat/X.flowers + (1|site/plot.treat), ziformula=~1, family="nbinom2",data = pot) #Potentilla model for negative binomial data with glmmadmb package
summary(pot.treat.glmmtmb) #summary of model output
## Family: nbinom2 ( log )
## Formula: total.seeds ~ treat/X.flowers + (1 | site/plot.treat)
## Zero inflation: ~1
## Data: pot
##
## AIC BIC logLik deviance df.resid
## 480.5 494.4 -232.2 464.5 34
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## plot.treat:site (Intercept) 3.090e-09 5.559e-05
## site (Intercept) 1.307e-02 1.143e-01
## Number of obs: 42, groups: plot.treat:site, 6; site, 4
##
## Overdispersion parameter for nbinom2 family (): 1.79
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.874253 0.573392 6.757 1.41e-11 ***
## treatopen-hand 0.636866 0.720335 0.884 0.3766
## treatopen:X.flowers 0.055490 0.032033 1.732 0.0832 .
## treatopen-hand:X.flowers 0.009497 0.029244 0.325 0.7454
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.5791 0.6073 -4.247 2.17e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
pot.counts.treat.simulation<-simulateResiduals(fittedModel = pot.treat.glmmtmb) #creating simulation for Potentilla mixed effects model with DHARMa package
plot(pot.counts.treat.simulation) #plotting simulation
testDispersion(pot.counts.treat.simulation) #checking for overdispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## ratioObsSim = 0.78442, p-value = 0.312
## alternative hypothesis: two.sided
testZeroInflation(pot.counts.treat.simulation) #checking for zero inflation
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 0.97784, p-value = 1
## alternative hypothesis: two.sided
The data are no longer zero-inflated, and we don’t have strong
evidence for pollen limitation in the Potentilla individuals.
Potentilla individuals are not pollen-limited.
Figures
Below are the plots for Potentilla individuals.
ggplot(pot.open, aes(x=pot.open$relative.position, y=pot.open$total.seeds, group=plot.treat)) +
labs(title="Effect of Potentilla Blooming Time on Total Developed Seeds", y="Total Developed Seeds",x="Days Relative to Population Peak Bloom") +
geom_point(aes(color=plot.treat),size=3, alpha = 0.8) +
scale_color_brewer(palette="Dark2") + labs(color="Plot Treatment") +
scale_color_manual(labels= c("control","accelerated"), values = c("control"="darkred","manip"="darkorange")) +
geom_smooth()
ggplot(pot.open, aes(x=pot.open$relative.position, y=pot.open$total.seeds, group=pot.open$early.late)) +
labs(y="Total Developed Seeds",x="Days Relative to Population Peak Bloom", title="Effect of Potentilla Blooming Time on Total Developed Seeds") +
theme_classic() +
geom_point(aes(color=early.late),size=3, alpha = 0.8) +
scale_color_brewer(palette="Dark2") + labs(color="Before or After Peak") +
scale_color_manual(labels= c("before peak","after peak"), values = c("early"="darkred","late"="darkorange")) +
geom_smooth(method = "glm", aes(color=early.late))
ggplot(pot.open, aes(x=pot.open$deviation, y=pot.open$total.seeds/pot.open$X.flowers, group= pot.open$early.late)) +
labs(y="Total Developed Seeds per Flower",x="Days Relative to Population Peak Bloom",title="Effect of Potentilla Deviation from Peak on Total Developed Seeds") +
theme_classic() +
scale_x_continuous(breaks=seq(0,24,2)) +
geom_point(size=3, alpha = 0.8, aes(color=early.late)) +
geom_smooth(method="glm",aes(color=early.late)) +
scale_color_brewer(palette="Dark2") +
labs(color="Before or After Peak") +
scale_color_manual(labels= c("before peak","after peak"), values = c("early"="darkred","late"="darkorange"))
ggplot(pot.open, aes(x=pot.open$number.conspecifics, y=pot.open$total.seeds/pot.open$X.flowers, group= pot.open$plot.treat)) +
labs(y="Total Developed Seeds per Flower",x="Conspecific Density", title="Effect of Potentilla Conspecific Density on Total Developed Seeds") +
theme_classic() +
geom_point(size=3, alpha = 0.8, aes(color=plot.treat)) +
geom_smooth(method="glm",aes(color=plot.treat)) +
scale_color_brewer(palette="Dark2") +
labs(color="Plot Treatment") +
scale_color_manual(labels= c("control","accelerated snowmelt"), values = c("control"="darkgreen","manip"="blue"))
Again we plotted the pollinator activity period.
ggplot(pot,aes(x=as.Date.factor(pot$midpoint), y=pot$total.seeds/pot$X.flowers, group= pot$treat)) +
labs(y="Total Developed Seeds per Flower",x="Date", title="Pollinator Activity Period and Total Developed Seeds in Potentilla") +
theme_classic() +
geom_point(size=3, alpha = 0.8, aes(color=treat)) +
geom_smooth(aes(color=treat)) +
scale_color_brewer(palette="Dark2") +
labs(color="Open or Hand-pollinated") +
scale_color_manual(labels= c("open","open-hand"), values = c("open"="darkred","open-hand"="darkorange")) +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
annotate("rect",xmin=as.Date.character("2019-07-05"), xmax=as.Date.character("2019-08-24"),ymin=0,ymax=Inf,alpha=0.1) +
scale_x_date(date_breaks = "3 days", date_labels = "%m/%d")
